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ABSTRACT 

The  geodesic  center  of  a  simple  polygon  is  a  point  inside  the 
polygon  which  minimizes  the  maximum  internal  distance  to  any 
point  in  the  polygon.  We  present  an  algorithm  which  calculates 
the  geodesic  center  of  a  simple  polygon  with  n  vertices  in  time 

1.   Introduction 

The  problem  addressed  in  this  paper,  to  locate  the  point  inside  a  simple 
polygon  P  whose  maximal  internal  distance  (by  a  route  inside  P)  from  any 
point  inside  P  is  minimal,  is  a  generalization  of  the  Euclidean  facility  location 
problem  which  asks  for  the  location  of  the  point  (facility)  which  is  least  far 
(in  the  Euclidean  metric)  from  the  furthest  of  a  finite  set  of  points  (the 
community  which  the  facility  is  to  serve)  [Me],  [Dy].  Indeed,  since  the 
furthest  point  from  a  point  inside  a  polygon  is  always  a  vertex  of  the 
polygon,  the  geodesic  center  of  the  convex  hull  of  the  community  to  be 
served  is  the  solution  to  the  standard  facility  location  problem.  We  can 
consider  the  problem  of  finding  the  geodesic  center  as  another  kind  of 
constrained  facility  location  problem  where  we  want  e.g.  to  locate  an 
emergency  service  on  a  polygonal  island  or  a  nurses  station  on  a  polygonal 
hospital  floor.  See  Fig.  1  for  an  illustration  of  the  geodesic  center  problem. 
The  standard  Euclidean  facility  location  problem  can  be  solved  in  time  0{n) 
([Me],  [Dy]),  but  its  extension  to  the  problem  of  finding  the  geodesic  center 
of  a  simple  polygon  appears  to  be  more  difficult. 
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The  problem  of  computing  the  geodesic  center  of  a  simple  polygon  has 
been  considered  by  T.  Asano  and  G.  Toussaint  [AT].  They  show  that  the 
geodesic  center  is  unique  and  present  an  algorithm  to  compute  it  in  time 
0{nHog  n),  where  n  is  the  number  of  vertices  in  the  given  polygon.  The 
main  idea  of  their  algorithm  is  to  construct  the  geodesic  furthest-point 
Voronoi  diagram  of  the  vertices  of  the  polygon  and  then  to  locate  the 
geodesic  center  at  either  a  vertex  of  the  Voronoi  diagram  or  at  the  midpoint 
of  a  geodesic  diameter  (i  e.  s  shortest  path  inside  the  polygon  joining  two 
vertices  which  has  maximal  length  over  all  choices  of  pairs  of  vertices).  We 
will  also  use  the  term  "geodesic  diameter"  to  denote  the  length  of  that  path. 
There  have  been  many  algorithms  to  find  the  geodesic  diameter  of  a  simple 
polygon.  The  best  result  at  the  present  time  is  an  0(n  log  n)  time  and  0(n) 
space  algorithm  due  to  S.  Suri  [Su]. 

A  related  problem  is  to  compute  the  link  diameter  and  the  link  center  of 
a  simple  polygon,  where  the  link  distance  between  two  points  is  the  minimum 
number  of  edges  in  a  polygonaa  path  joining  them  inside  the  polygon  and 
where  the  link  center  and  diametei  are  defined  in  an  analogous  manner  to  the 
definition  of  geodesic  center  drid  diameter.  In  this  case  the  link  center  is  no 
longer  unique  but  consists  of  a  polygon  wliich  may  be  as  large  as  the  entire 
given  polygon.  Suri  [Su2]  has  an  0(n  iog^/i)  time  and  0{n)  space  algorithm 
which  computes  the  hnk  diameter  of  a  simple  polygon  and  [Lea]  presents  an 
O(n^)  algorithm  for  computing  the  link  center  of  a  simple  polygon.  Suri  and 
El-Gindy  (private  communication)  also  report  similarly  efficient  algorithms 
for  computing  the  link  center. 

Our  algorithm  proceeds  as  follows.  We  start  with  a  triangulation  of  the 
polygon  P,  then  perform  something  like  a  binary  search  through  the 
diagonals  of  the  triangulation,  determining  at  each  tested  diagonal,  via  the 
algorithm  RELCEN  to  be  describe  in  Section  3  below,  on  which  side  of  that 
diagonal  the  geodesic  center  lies.  In  this  way  we  locate  a  triangle  which 
contains  the  geodesic  center.  However,  as  we  move  around  in  this  triangle 
the  combinatorial  structure  of  the  shortest  path  to  a  given  vertex  can  change. 
We  next  use  a  modification  of  Megiddo's  method  for  doing  linear 
programming  in  linear  time  [Me]  to  find  a  polygon  R  within  this  triangle 
where  the  structure  of  the  path  from  each  point  x  ^  R  to  each  vertex  v  of  P 
remains  constant,  and  consists  of  the  straight  segment  from  x  to  some, 
possibly  different,  vertex  followed  by  a  fixed  path  of  known  length  to  v.  In 
other  words,  in  our  final  subregion  R,  the  internal  distance  from  a  point  x  to 
the  j-th  vertex  of  P  has  a  fixed  analytic  expression  of  the  form  [xw,|+c,,  for 
1  =  1,2,  .  .  .  ,n.  The  problem  that  remains  is  equivalent  to  that  of  finding  a 
smallest  circle  that  contains  a  given  collection  of  n  circles  (where  the  z-th 
circle  is  centered  at  u,  and  has  radius  c,).   This  final  task  is  accomplished  by 
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another  technique  of  Megiddo  [Me2]  for  deriving  efficient  sequential 
optimization  algorithms  from  parallel  algorithms. 

The  total  time  complexity  of  our  algorithm  is  0{n  log^n).  The  final  stage 
(that  of  finding  the  smallest  spanning  circle  of  circles),  which  may  be  of 
independent  interest  in  location  theory,  runs  in  time  0{n  log  n  log  log  n). 

The  paper  is  organized  as  follows.  In  section  2  we  present  the  geometric 
preliminaries  needed  to  justify  the  algorithm  RELCEN  which  is  presented  in 
section  3.  In  section  4  we  present  the  complete  algorithm,  and  concluding 
remarks  are  given  in  Section  5. 

i 
2.  Geometric  Preliminaries 

Definitions.  Let  P  be  a  simple  polygon  with  n  vertices  (more  precisely  P  is  a 
closed  simply-connected  planar  region  whose  boundary  is  a  given  simple 
polygon).  For  x,y  in  P,  g(x,y)  is  the  shortest  path  inside  P  between  x  and  y, 
d{x,y)  is  the  length  of  this  path.  Furthermore  let  xy  denote  the  straight 
segment  with  endpoints  x  and  y  and  \xy\  its  length.  Finally,  we  let  u{x,y) 
denote  the  unit  vector  in  the  direction  that  the  path  gix,y)  starts  at  from  x, 
and  define  the  angle  Ucyz  as  the  smallest  of  the  two  angles  between  u(y,x) 
and  u(y,z).  A  subset  C  of  P  is  called  P -convex  if  whenever  x,y  are  in  C, 
g(x,y)  is  a  subset  of  C.  It  immediately  follows  that  the  intersection  of  P- 
convex  sets  is  P-convex,  and  that  P-convex  sets  are  connected. 

As  shown  in  Lee  and  Preparata  [LP],  g{x,y)  is  a  polygonal  path  whose 
corners  are  vertices  of  P.  For  a  fixed  point  x  in  P,  the  union  of  g{x,v)  over 
all  vertices  v  of  P  is  a  planar  tree  Q{x)  (rooted  at  x),  which  we  call  the 
shortest  path  tree  of  P  (with  respect  to  x).  This  tree  has  n  nodes,  namely  the 
vertices  of  P,  and  its  edges  are  straight  segments  connecting  these  nodes.  It 
has  been  shown  by  Guibas  et  al.  [Gea]  that  this  tree  can  be  computed  in 
linear  time. 

Given  three  points  a,b,c  in  P,  Consider  the  geodesic  paths  g{a,b), 
g{b,c)  and  g(c,a).  There  exist  points  a' ,  b'  and  c'  such  that  the  paths  g{a,b) 
and  g{a,c)  intersect  in  the  path  g(a,a'),  the  paths  g{b,c)  and  g{b,a)  intersect 
in  g{b,b')  and  finally  the  paths  gic,a)  and  gic,b)  intersect  in  g{c,c'). 
Moreover  we  see  that  the  geodesic  triangle  Aa'b'c'  has  only  reflex  angles 
along  its  boundary  with  the  exception  of  the  angles  at  a' ,  b'  and  c'  (see  Fig. 
2). 

Lemma  L  As  a;  varies  openly  along  g(b,c),  d(a^)  is  a  unimodal  function 
which  has  its  maximum  at  fc  or  at  c,  i.e.  d(a,x)  <  max  {d(a,b),  d{a,c)}. 

Proof.  Without  loss  of  generality  we  assume  that  ^.abc  is  a  geodesic  triangle, 
i.e.  g(a,b),  g{b,c)  and  g{a,c)  are  openly  disjoint,  and  that  g{a,x)  and  g{a,b) 
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are  also  openly  disjoint.  Let  q  be  the  last  point  that  g{a^)  and  g{a,c)  have  in 
common.  Let  x'  and  x' '  be  the  intersection  of  qx  with  ab  and  be  respectively 
(see  Fig.  2).  Since  d{a,q)  <\ax'\+  \x:'q\,  we  have  d{a^)  <  |ax'|  +  \k'x"\. 
If  [x'a;"|  <  [x'^l,  then  d{a^)  <  \ab\  <  d(a,b).  Otherwise  the  angle  Lqx"c  is 
certainly  >  tt/I,  whence  \qx"\  <  \qc\.   Hence 

dia^)  <  d(a,q)  +  \qx"\<  d{a,q)  ^  \qc\<  d{a,q)  +  d(q,c)  =  d(a,c)  . 

Q.E.D. 

Corollary  1.  The  closed  ball  of  radius  r  at  x,  B(r^)  =  {y  ^  P  \  d{x,y)  <  r}, 
is  P-convex. 

The  proof  of  the  lemma  also  implies 

Corollary  2.  With  the  notation  of  the  lemma,  if  the  angle  Lb'a'c'  at  a'  is 
greater  than  or  equal  to  7t/2  then  d{b,c)  >  d(a,b),  d{a,c). 

It  is  obvious  that  d{a,x)  is  maxiraired  for  x  a  vertex  of  P. 

Let  F{x)  =  max  d{x,v),  where  ihe  maximum  is  taken  over  the  vertices  v 
of  P.  Qearly  F  is  continuous  and  thercfoi'e  has  a  minimum  value  called  the 
geodesic  radius  of  P  and  denoted  gr(P).  The  set  of  points 

geocen{P)  =  {x  ?  P  \F{x)  =  gr{P)} 

is  the  intersection  of  the  closed  bells  of  radius  gr{P)  centered  at  the  vertices 
of  P.  If  geocen{P)  contained  two  distinct  points  x,y  it  would  have  to  contain 
g{x,y)  and  the  geodesic  distance  from  a  vertex  would  have  to  be  constant 
along  each  segment  of  g{x,y),  which  is  clearly  impossible  by  Lemma  1. 
Hence  geocen{P)  is  a  singleton  called  the  geodesic  center  of  P  (cf.  also 
[AT]). 

Let  c  be  a  point  of  P  and  let  the  furthest  vertices  from  c  be 
{vi,V2,  .  .  .  ,v^}  =  V(c)  and  suppose  that  they  all  have  geodesic  distance  r 
from  c.  For  any  v  in  V(c),  the  line  orthogonal  to  u(c,v)  determines  a 
diagonal  D  of  P  which  separates  P  into  two  components  (the  "sides"  of  D). 
Suppose  that  j:  is  a  point  of  P  lying  on  the  side  of  D  which  does  not  contain 
v;  then  g{x,v)  must  cross  D,  say  at  x'  and  since  the  angle  Ljc'cv  is  tt/Z  we  see, 
by  Corollary  2,  that  d{x,v)  >  d{x'  ,v)  >  d{c,v)  =  F{c),  whence  x  is  not 
geocen{P)  (see  Fig.  3).   As  a  consequence  we  have 

Lemma  2.  If  the  set  {u(c,v)  [  v  6  V(c)}  does  not  lie  in  a  an  open  half-plane 
through  c  then  c  =  geocen(P). 

Remark.  If  two  directions  u{c,v{)  and  u(c,V2)  are  opposite,  then  again  c  is 
the  geodesic  center  and  g(vi,v2)  is  a  geodesic  diameter  (whose  midpoint  is  c). 

Definition.  Given  a  cut  of  P,  i.e.  a  segment  C  which  is  openly  contained  in 
P  and  has  endpoints  on  the  boundary  of  P,  thus  separating  P  into  two 
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components,  the  relative  center  of  P  on  C  is  the  unique  point  c  of  C  which 
minimizes  F{x)  for  x  in  C  and  is  denoted  RELCEN{C)  =  c. 

Lemma  3.  Let  c  =  RELCEN{C),  and  suppose  that  all  directions  u(c,v)  for  v 
in  V(c)  lie  in  an  open  half -plane  bounded  by  the  line  H  through  c.  Let  u*  be 
a  unit  vector  orthogonal  to  H  and  in  this  half-plane.  Then  geocen(P)  lies  on 
the  side  of  C  to  which  u*  points  (see  Fig.  4;  the  proof  will  imply  that  u* 
cannot  point  along  C). 

Proof.  Let  d{c,v)  =  r  for  each  v  in  V(c),  while  the  next  furthest  vertex  from 
c  has  geodesic  distance  r'  <  r  from  c.  Note  that  for  c'  in  B{r—r',c)  the 
distance  from  c'  to  any  vertex  not  in  V(c)  'is  at  most  r.  Choose  a  positive 
8  <  min  {r— r',p},  where  p  is  the  length  of  the  smallest  initial  segment  of 
g{c,v)  for  V  in  V(c),  sufficiently  small  so  that  the  boundary  oi  B  =  B(8,c) 
consists  of  a  single  circular  arc  and  possibly  a  portion  of  at  most  two 
segments  of  the  boundary  of  P  (the  latter  might  occur  in  case  c  is  on  the 
boundary  of  P).  For  each  v  6  V{c)  let  /z(v)  be  the  peri^endicular  distance 
from  H  to  py,  where  p^  is  the  intersection  of  g{c,v)  and  the  boundary  of 
5(8, c),  and  let  h  =     min    h(v).    The  point  c  +  hu*  =  p   in  5(8, c)  -  C 

V  €  V(c) 

satisfies  d(p  ,v)  <  r  for  all  vertices  v  of  P,  since  for  v  in  V(c), 

d(p,v)  <  d{p,py)  +  d(py,v)  <  d{c,py)  -H  d(py,v)  =  r  , 

and  we  have  already  seen  that  d(p,v)  <  r  for  the  other  vertices  of  P.  Hence 
p  lies  in  the  intersection  of  the  closed  balls  of  some  radius  r"  <  r  at  each  of 
the  vertices  of  P.  Since  the  geodesic  center  of  P  must  also  lie  in  this 
intersection  which  is  connected  (since  P-convex)  and  disjoint  from  C  (since, 
by  definition  of  RELCEN{C),  F  >  r  >  r"  on  C),  we  see  that  the  geodesic 
center  of  P  lies  on  the  same  side  of  C  as  p.  Q.E.D. 

3.   The  algorithm  RELCEN 

Let  C  be  a  cut  of  P,  i.e.  C  is  openly  contained  in  P  and  has  endpoints  a 
and  b  on  the  boundary  of  P.  We  construct,  in  linear  time  by  [Gea],  the  tree 
Q{a)  of  shortest  paths  from  a  to  the  vertices  of  P  and  the  tree  Q{b)  of 
shortest  paths  from  b  to  the  vertices  of  P.  For  each  edge  e  of  either  of  these 
trees  we  find  the  intersection  (e,C)  of  the  line  determined  by  e  and  the  line 
determined  by  C.  If  (e,C)  lies  in  C  and  the  segment  from  (e,C)  to  e  lies  in 
P  (the  algorithms  in  [Gea]  provide  this  information  for  all  edges  e  in 
Q(a)  U  Q{b)  in  overall  linear  time),  then  we  mark  the  point  (e,C)  on  C  (see 
Fig.  5).  These  0(n)  points  are  critical  points  where  the  structure  of  the 
shortest  path  from  a  point  x  on  C  to  some  vertex  v  of  P  and  the  analytic 
expression  for  d{x,v)  change;  i.e.  for  j:  in  a  subinterval  determined  by  these 
0(n)   points,   the  analytic  form   of  d(x,v)   is   uniquely   determined  by  v. 
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Namely,  d{x,v)  =  [xu^j  +  Cy,  where  Uy  is  a  vertex  of  P  determined  by  v.  Our 
object  is  next  to  locate  the  subinterval  of  C  in  which  RELCEN{C)  must  lie. 
The  functions  d{x,v)  for  v  a  vertex  of  P  are  unimodal  on  C,  hence 
F{x)  =  max  d{x,\>)  is  also  unimodal.  Given  a  point  c  of  C,  we  can  determine 

V 

the  side  of  c  on  which  RELCEN{C)  is  located  by  again  calculating  (in  0{n) 
time)  the  tree  Q{c)  of  shortest  paths  from  c  to  the  vertices  of  P  and  then  by 
finding  the  furthest  vertex  v  from  c.  Then  F{c)  =  d{c,v),  and  the  side  of  c 
containing  RELCEN{C)  (which  is  the  point  where  F{x)  is  minimized  on  C)  is 
determined  by  the  sign  of  ttie  derivative  of  d{x,v)  at  or  =  c  (a  similar  rule 
japplies  when  F{c)  is  simultaneously  attained  by  more  than  one  vertex  v). 
Hence,  performing  a  binary  search  over  the  0{n)  critical  points  (e,C)  on  C, 
we  can  find  in  time  0{n\ogn)  a  subinterval  [a',fc']  of  C  in  which 
RELCEN{C)  lies,  and  over  which  the  analytic  form  of  d{x,v),  for  every 
vertex  v  of  P,  is  determined  solely  by  v. 

Without  4oss  of  generality  we  may  assume  that  we  are  now  on  the 
interval  C  =  [0,1],  over  which  we  are  given  n  functions  di  of  the  form 

d^x)  =   [{x-xiy  i-  yf]       +  c,  ,       x^  [0,1]  , 

for  i  =  I,  .  .  .  ,n,  and  we  want  to  locate  the  unique  minimum  of  the  upper 
envelope  F{x)  of  these  functions. 

This  final  step  of  the  algorithm  is  substantially  simpler  than  the 
preceding  steps,  and  can  be  implemented  in  time  linear  in  the  number  of 
functions  di;  note  also  that  the  goal  of  this  step  is  stated  in  a  manner 
independent  of  the  structure  of  P.  The  method  we  use  is  a  variant  of  the 
linear  time  technique  of  Megiddo  for  linear  programming  in  R^  [Me],  and 
runs  as  follows. 

We  group  the  functions  di  in  at  most  n/2  disjoint  pairs,  so  that  the  two 
functions  in  each  pair  correspond  to  two  vertices  of  P  that  lie  on  the  same 
side  of  C.  Let  di,  dj  be  such  a  pair.  We  claim  that  these  functions  have  at 
most  one  intersection  point  on  C.  Indeed,  suppose  to  the  contrary  that 
di{x)  =  dj{x)  and  diiy)  =  dj(y)  for  two  distinct  points  x,y  i  C.  Let 
v/  =  {xi,yi),  vj  =  (xj,yj).  Then  we  have 

[xvj   +    Ci    =    \xVj\   +    Cj 

\yv,\  +  c,  =  \yvj\  +  Cj  , 

or 

\xvj\  +  \yvi\  =  [xv^I  +  \yvj\  . 

Since  v,,  Vj  are  vertices  of  P  lying  on  the  same  side  of  C,  and  the  four 
segments  xvi,  xvj,  yvj,  yvj  are  all  contained  in  P,  it  follows  that  either  xvj  and 


yvi  intersect,  or  xv,  and  yvj  intersect,  or  the  segments  xv,,  xvj  overlap,  or  the 
segments  yvt,  yvj  overlap.  In  the  first  case  we  have,  using  the  triangle 
inequality  (cf .  Fig.  6) 

\kvj\  +  \yvt\  >  \xvi\  +  \yvj\  , 

which  contradicts  the  previous  equality.  Similar  contradictions  are  obtained  in 
each  of  the  remaining  cases,  thus  establishing  our  claim. 

We  collect  these  0{n)  critical  points,  and  calculate  their  median  xq  in 
0(n)  time.  We  then  find,  as  above,  which  function(s)  attain  the  maximum  F 
at  Xq,  and  either  determine  that  xq  is  the  relative  minimum  of  F  on  C,  or  else 
find  which  side  of  xq  contains  this  relative  minimum  c.  Thus,  for  each  of  the 
n/4  pairs  of  functions  whose  intersection  lies  on  the  other  side  of  xq,  we 
know  which  of  the  two  functions  dominates  the  other  at  c ,  and  can  therefore 
discard  the  smaller  function  from  further  consideration.  Thus,  in  0(n)  time 
we  discard  a  quarter  of  the  functions  d,,  and  then  repeat  the  same  step  for  the 
remaining  3n/4  functions.  After  at  most  0{\og  n)  iterations,  the  number  of 
remaining  functions  becomes  sufficiently  small  to  allow  explicit  calculation  of 
c  in  constant  time.  Thus  the  entire  step  runs  in  linear  time  (see  [Me]  for  more 
detail). 

Let  c  be  the  relative  center  of  P  on  C.  Again  we  construct  in  0{n)  time 
the  tree  of  shortest  paths  Q{c)  from  c  to  the  vertices  of  P.  Lemmas  2  and  3 
allow  us  to  determine  either  that  c  is  the  geodesic  center  of  P,  or  else  on 
which  side  of  C  that  geodesic  center  must  lie. 

It  is  clear  that  the  complexity  of  RELCEN  is  0(n  log  n).  The  algorithm 
RELCEN  is  used  as  a  subroutine  in  the  complete  algorithm  GEOCEN  given 
below,  with  the  convention  that  if  the  output  produced  by  RELCEN  is  the 
geodesic  center  of  P  then  the  entire  algorithm  halts  and  outputs  this  center. 

4.  The  algorithm  GEOCEN. 

We  are  now  in  position  to  describe  the  complete  procedure  for  the 
calculation  of  geocen{P).  It  consists  of  the  following  steps. 

I.  Triangulate  P  (in  time  0{n  log  n)  using  the  algorithm  in  [GJPT],  or  more 
efficiently  using  the  algorithm  of  Tarjan  and  Van  Wyk  [TV]),  and  place  the 
internal  edges  of  the  triangulation  in  a  balanced  binary  tree  T  [Ch]  such  that 
for  each  subtree  T  of  7,  the  edges  stored  at  the  nodes  of  T  bound  a 
collection  of  triangles  that  cover  a  connected  simple  subpolygon  of  P,  and 
such  that  the  number  of  edges  in  each  of  the  two  subtrees  rooted  at  the  left 
and  the  right  children  of  a  node  e  of  T  are  both  at  least  one  third  of  the 
number  of  edges  in  the  subtree  rooted  at  e.  Such  a  balanced  decomposition 
tree  of  P  can  be  constructed  in  linear  time  [Gea],  or,  by  an  alternative 
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technique,  in  0{n  log  n)  time  [Ch]. 

We  next  search  the  tree  T  to  locate  the  triangle  Aabc  of  the  triangulation  of  P 
which  contains  the  geodesic  center  of  P.  At  each  node  of  T  along  the  search 
path,  we  apply  RELCEN  to  the  corresponding  diagonal  e  of  P  to  determine 
which  side  of  e  contains  geocen{P),  and  then  continue  the  search  at  the 
appropriate  child  of  that  node.  Hence  in  0(log  n)  applications  of  RELCEN, 
each  of  which  takes  0(n  log  n)  time,  we  obtain  the  desired  triangle.  Thus  the 
total  time  required  by  this  step  is  0{n  log^/i). 

n.  Next  we  want  to  locate  a  region  within  the  triangle  Aabc  in  which,  for 
each  vertex  v  of  F,  the  structure  of  the  path  g(x,v)  and  the  analytic  form  of 
the  corresponding  function  d{x,v)  sre  both  uniquely  determined  by  the  vertex 

V. 

To  this  end  we  construct  the  trees  Q{a),  Q{b),  Q(c)  of  shortest  paths  from 
a,  b  and  c  respectively,  and  record  ihe  k  =  0(c)  "tangents",  i.e.  extensions 
of  shortest  path  edges  toward  the  triangle  Aabc  that  actually  enter  that 
triangle  (see  Fig.  7).  (As  a  matter  of  fact,  we  only  need  to  calculate  the 
portion  of  the  tree  Q(a)  ccrxsisting  of  shortest  paths  that  do  not  cross  the 
opposite  side  be  of  our  triangle,  and  sirnil&jly  for  Q(b),  Q{c).)  Within  any 
region  determined  by  these  tangents  ihe  structure  of  g(x,v)  and  the  analytic 
form  of  d(x,v)  are  both  uniquely  determined  by  the  vertex  v.  Our  task  is  thus 
to  determine  the  side  of  each  one  of  these  tangents  which  contains  geocen(P). 
To  this  end,  we  use  and  adapt  a  technique  used  by  Megiddo  [Me]  in  his 
linear-time  algorithm  for  linear  programming  in  R-'.  Specifically,  we  find  the 
median  slope  of  these  tangents  and  let  T-^,  .  .  .  ,T^f2  have  positive  slope  and 
Tkn+\>  •  •  •  'Tk  have  negative  slope  with  respect  to  a  new  orientation  of  the 
axes  where  the  x-axis  is  in  the  direction  of  the  median  slope.  Form  the  k/2 
intersections  of  the  lines  T,  and  7^/2+ 1  for  i  =  1,  .  .  .  .k/2.  Next  find  the 
median  x^  of  the  j: -coordinates  of  these  k/2  intersection  points.  Take  the  cut 
Cj  determined  by  the  line  x  =  j^i  in  the  triangle  Aabc  and  apply  RELCEN  to 
Cjc  to  determine  which  side  of  C^  contains  the  geodesic  center.  For  the  k/4 
intersection  points  lying  on  the  other  side  of  C^  find  their  median  y- 
coordinate  >'i.  Take  the  cut  Cy  determined  by  the  line  y  =  ji  in  this  triangle 
and  apply  RELCEN  to  Cy  to  determine  which  side  of  Cy  contains  the 
geodesic  center.  Now  for  each  of  the  it/8  pairs  of  tangents  whose  intersection 
lies  in  the  quadrant  determined  by  C^  and  Cy  and  opposite  to  the  quadrant 
containing  the  center,  we  can  discard  one  of  the  two  tangents  (since  the  side 
of  it  containing  geocen{P)  is  now  known).  In  addition,  we  add  the  two  half- 
planes  bounded  by  the  cuts  C^  and  Cy  and  containing  geocen{P)  to  a 
collection  of  half-planes  that  the  algorithm  constructs.  Thus  in  C>(log  n) 
applications  of  RELCEN  we  can  find  the  side  of  each  tangent  containing 
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geocen{P),  and  at  the  same  time  obtain  0(log  n)  half -planes  whose 
intersection  is  known  to  contain  geocen{P).  We  then  find  an  arbitrary  point 
w  lying  in  the  intersection  of  these  half-planes,  a  task  that  can  be 
accomplished  in  0(log  n)  time,  using  Megiddo's  technique  for  two- 
dimensional  linear  programming  [Me],  and  then  calculate  the  shortest-path 
tree  Q{w)  in  0{n)  time.  By  the  preceding  discussion,  the  structure  of  Q{w) 
and  of  Q{geocen{P))  are  the  same,  so  that  the  analytic  form  of  each  function 
d{x,v)  throughout  the  intersection  of  those  half-planes  can  be  readily 
determined. 

Remark.  In  the  preceding  arguments  we  have  assumed  that  geocen(P)  is  an 
interior  point  of  P.  It  is  easy  to  check,  using  Lemma  2,  that  if  geocen{P)  lies 
in  the  relative  interior  of  an  edge  of  P,  the  procedure  sketched  above  will 
locate  a  region  containing  the  center,  and  will  properly  obtain  the  fixed 
'  • :'  tic  form  of  the  functions  d(x,v)  for  x  in  that  region.  Finally,  if 
geycen{P)  is  a  vertex  of  P  (which  is  to  say,  geocen{P)  is  a,  b,  or  c)  then  it 
must  be  a  reflex  vertex,  and  so  be  incident  to  a  side  of  Aabc  which  is  also  a 
cut  of  P.  But  then  a  preceding  application  of  RELCEN  to  that  cut  would 
have  found  geocen{P)  to  lie  on  that  cut,  and  consequently  would  have 
terminated  our  algorithm.  Hence  in  step  EI,  and  in  the  suc-ceeding  step  lU,  we 
can  assume  that  geocen{P)  is  not  a  vertex  of  P. 

in.  The  problem  is  now  reduced  to  finding  the  minmax  of  a  family  of 
functions  di{x)  =  d{x,Vi)  =  ixv^(/)|  +  C/,  where  C/  =  ^(Vit(/),v,),  for 
i  =  1,  .  .  .  ,n.  Note  that  the  preceding  remark  and  discussion  allow  us  to 
ignore  P  completely  and  regard  each  of  the  functions  di{x)  as  being  defined 
over  the  entire  plane,  because  geocen{P)  must  be  the  unique  global  minimum 
of  F{x)  =  max  di{x). 

This  final  task  can  be  accomplished  by  a  technique  due  to  Megiddo  [Me2] 
which  uses  parallel  algorithms  in  the  design  of  efficient  sequential 
optimization  algorithms.  To  apply  this  technique  we  first  consider  a 
restricted  version  of  the  problem,  in  which  we  want  to  determine  the  relative 
minimum  F{c)  of  F{x)  on  a  given  Hne  C,  where  c=RELCEN{C),  and  also 
determine  on  which  side  of  C  the  global  minimum  of  F  lies.  This  problem 
can  in  fact  be  solved  by  a  similar  procedure  to  the  last  step  of  the  RELCEN 
algorithm,  which  runs  in  time  proportional  to  the  number  of  functions  (and  is 
independent  of  the  number  of  vertices  of  P!).  For  exposition  sake  we  sketch 
here  the  variant  of  the  last  step  of  RELCEN,  denoted  here  by  RELCEN*, 
needed  in  this  case.  It  consists  of  0(log  n)  substages,  in  each  of  which  a  fixed 
fraction  of  the  functions  di(x)  still  present  is  being  discarded  (because  J,(c)  is 
determined  to  be  less  than  F(c)),  until  we  are  left  with  only  0{1)  functions, 
whose  relative  minmax  on  C  (as  well  as  the  side  of  C  containing  the  global 
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minimum  of  F(x))  can  be  found  in  constant  time.  Each  such  substage 
combines  the  remaining  k  functions  in  disjoint  pairs,  determines  for  each  pair 
di,  dj  the  (at  most  two)  points  on  C  at  which  di  =  dj,  thus  obtaining  at  most 
k  critical  points  on  C.  We  then  find  the  median  xq  of  these  points,  and 
calculate  in  0{k)  time  the  side  of  xq  containing  the  relative  minmax  c  on  C, 
by  calculating  F(xo)  =  max  di{xo)  and  evaluating  the  derivative  of  F  along  C 

at  Xq.  Next  we  find  the  median  xi  of  the  critical  points  lying  on  the  side  of  xq 
containing  c,  and  again  find  the  side  of  x-^  containing  c.  Thus  we  can 
determine,  for  at  least  U4  of  the  pairs  of  functions,  which  of  the  two 
functions  in  the  pair  dominates  the  other  at  c,  and  thus  discard  the  smaller 
one.  The  remaining  3/:/4  functions  are  then  passed  to  the  next  substage.  We 
continue  in  this  manner  imtil  the  number  of  remaining  fimctions  is 
sufficiently  small  to  allow  direct  calci.ilation  of  c  in  0(1)  time.  (The  idea  of 
performing  two  median  calculations  at  e-ach  step  of  RELCEN*  is  taken  from 
[Ze].)  Qearly  the  time  complexity  of  RE-LCEN*  is  0{n). 

To  find  the  desired  glob.il  minmax  z*  =  {x*,y*),  we  simulate  the 
execution  of  RELCEN*  on  the  line  y~y*,  without  knowing  explicitly  the 
value  of  y*.   Denote  this  line  by  C. 

The  first  step  of  REI.CErs'"  is  to  group  the  n  functions  di{x)  in  n/2 
disjoint  pairs.  Each  pair  d^,  dj  determines  a  simple  unbounded  curve  7,^ 
which  is  either  a  branch  of  a  hyperbola  or  a  straight  line,  on  which  di  =  dj. 
Each  ytj  intersects  C  in  at  most  two  points,  unless  ■y,j  coincides  with  C;  we 
will  assume  that  this  latter  case  does  not  arise;  it  can  be  handled  by  minor 
modifications  of  the  algorithm  described  below.  Before  proceeding,  we  have 
to  determine  exactly  at  how  many  points  each  7^  intersects  (the  still 
unknown)  C.  This  is  done  as  follows.  If  7,^  does  not  have  a  point  of 
horizontal  tangency,  then  it  must  intersect  C  in  precisely  one  point. 
Otherwise  the  number  of  intersection  points  is  either  0  or  2,  depending  on 
the  position  of  C  relative  to  the  horizontal  tangent  of  7^.  We  thus  collect  the 
0{n)  y -coordinates  of  these  tangents,  and  locate  y*  among  them  by  a  binary 
search,  each  step  of  which  is  performed  by  applying  RELCEN*  to  the 
appropriate  horizontal  line.  Thus  in  total  0(n  log  n)  time  we  can  determine 
how  many  intersections  C  has  with  all  the  hyperbolic  branches  7/y. 

The  next  step  of  RELCEN*  is  to  calculate  the  median  of  those  (at  most 
n)  critical  points.  To  do  this  we  execute  a  fast  parallel  median-finding 
algorithm.  The  asymptotically  best  algorithm  of  this  kind  is  that  by  Ajtai, 
Komlos,  Steiger  and  Szemeredi  [AKSS]  which  runs  in  C>(log  log  n)  parallel 
time  and  uses  0(n)  processors. 

Continuing  our  algorithm,  we  consider  one  parallel  step  of  the  median- 
finding  algorithm.    It  performs  0{n)  comparisons  between  pairs  of  critical 
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points.  For  each  such  comparison,  the  corresponding  pair  of  hyperbolic 
branches  intersects  in  the  plane  in  at  most  four  points,  and  the  outcome  of 
the  comparison  depends  on  the  position  of  y*  relative  to  the  >'-levels  of  these 
points.  We  sort  the  collection  of  ^-coordinates  of  these  0(n)  intersection 
points,  and  perform  as  above  a  binary  search  to  locate  y*  among  them.  It 
thus  follows  that  in  0{n  log  n)  time  we  can  resolve  each  of  the  comparisons 
in  the  present  parallel  step  of  the  [AKSS]  algorithm,  so  that  the  entire 
sequential  simulation  of  this  algorithm  (at  the  still  unknown  level  y*)  can  be 
done  in  time  0(n  log  n  log  log  n). 

We  have  thus  foipd  the  first  median  xq  calculated  in  RELCEN*  (actually 
we  do  not  have  an  explicit  value  of  xq,  but  only  the  hyperbolic  branch  -y^  on 
which  the  (unknown)  point  {xQ,y*)  lies).  The  next  step  of  RELCEN*  is  to 
calculate  the  maximimi  of  the  n  values  dt{xQ,y*),  from  which  the  side  of  xq 
containing  the  relative  minmax  of  F  on  C  can  be  easily  determined.  To  this 
end  we  run  a  second  fast  parallel  algorithm,  this  time  for  finding  the 
maximum  of  n  numbers;  we  can  use  e.g.  the  algorithm  of  Valiant  p/a]  which 
also  runs  in  0(log  log  n)  parallel  time,  and  uses  0{n)  processors.  At  each 
parallel  step  of  this  algorithm  we  perform  0{n)  comparisons  of  the  form 
^tC^oJ*)  vs.  di{xQ,y*),  subject  to  the  constraint  that  {xQ,y*)  lies  on  the 
hyperbola  7^^,  calculated  by  the  preceding  median-finding  step.  It  follows  that 
each  such  comparison  determines  at  most  four  critical  j-levels  in  which  the 
result  of  the  comparison  can  change.  Collecting  all  those  0{n)  critical  y- 
levels,  and  performing  a  binary  search  over  them  as  before,  we  can  complete 
the  current  parallel  step  of  the  maximum  finding  algorithm  in  0{n  log  n) 
time,  and  thus  the  sequential  simulation  of  the  entire  algorithm  (at  the 
unknown  y*)  can  be  done  in  0(n  log  n  log  log  n)  time. 

At  this  point  we  have  to  calculate  the  sign  of  the  partial  A:-derivative  of  F 
at  {xQ,y*).  It  is  easily  checked  that  this  sign  can  change  in  at  most  two  critical 
J -levels,  so  again  we  can  determine  this  sign  by  at  most  two  applications  of 
RELCEN*  to  the  horizontal  lines  at  those  j -levels.  Next  we  have  to 
determine  which  of  the  critical  points  on  C  lie  on  the  appropriate  side  of  xq. 
This  is  achieved  by  comparing  each  critical  point  with  xq.  Again,  each  such 
comparison  determines  at  most  four  critical  j-levels  at  which  its  outcome  can 
change,  and  a  binary  search  through  those  0(n)  j-levels,  using  RELCEN*  as 
above,  will  resolve  all  necessary  comparisons  in  0(n  log  n)  time. 

The  next  steps  of  the  simulated  RELCEN*  consist  of  a  second  median 
calculation,  followed  by  a  second  maximum  calculation  and  a  second 
classification  of  critical  points  on  C,  all  of  which  can  be  simulated  as 
described  above. 
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Thus  in  time  0{n  log  n  log  log  /j)  we  can  complete  the  present  stage  of 
the  simulated  RELCEN*,  thereby  eliminating  n/4  of  the  functions  di.  It 
therefore  easily  follows  that  the  entire  simulation  of  RELCEN*  at  y  =  y* 
can  also  be  accomplished  in  time  0{n  log  n  log  log  n).  When  this  simulation 
terminates,  we  are  left  with  0(1)  functions  di,  whose  global  minmax  can  now 
be  determined  explicitly  in  constant  time.  Thus  this  final  stage  of  the 
GEOCEN  algorithm  can  be  accomplished  in  0{n  log  n  log  log  n)  time. 

Remark:  As  noted  in  the  introduction,  the  problem  solved  in  this  last  phase 
of  our  algorithm  can  be  restated  as  follows:  Given  n  circles  in  the  plane,  find 
the  smallest  circle  which  contains  all  of  them.  While  the  simpler  1-center 
problem  (in  which  instead  of  n  circles  we  are  given  n  points)  has  been  solved 
in  0{n)  time  ([Me],  [Dy]),  we  do  not  know  of  similarly  efficient  solutions  for 
the  case  of  circles.  Note  also  that  the  constants  involved  in  the  complexity 
bounds  of  the  algorithm  of  [AKSvS]  that  we  have  used  are  fairly  large.  We 
could  alternatively  use  simpler  parallel  techniques,  e.g.  the  algorithm  of  Cole 
and  Yap  [CY]  that  runs  in  parallel  time  0((log  log  n)^)  and  uses  0{n) 
processors,  or  even  a  simpler  sorting  scheme,  such  as  that  of  Preparata  [Pr] 
or  its  recent  improvement  by  Cole  [Co].  Use  of  these  simpler  schemes  will 
lead  to  more  efficient  algorithms  in  practice;  moreover,  our  entire  algorithm 
has  already  consumed  0(n  log^/i)  time  so  far,  and,  as  will  follow  from  the 
analysis  given  below,  use  of  any  of  these  parallel  schemes  will  still  yield  an 
0(n  log^n)  algorithm  for  the  final  phase  at  hand.  Nevertheless,  since  this 
final  phase  may  be  of  independent  interest  in  facility  location  theory,  we  have 
used  the  algorithm  of  [AKSS]  to  obtain  a  solution  to  that  minmax  problem 
with  the  best  asymptotic  complexity.) 

5.  Discussion 

We  have  presented  a  rather  involved  algorithm  for  computing  the 
geodesic  center  of  an  /i-sided  simple  polygon  in  time  0{n  log^n).  While  this 
constitutes  a  significant  improvement  over  the  previous  algorithm  in  [AT], 
we  have  no  non-linear  lower  bound  on  the  complexity  of  the  problem.  An 
obvious  open  problem  is  whether  the  geodesic  center  can  be  calculated  in 
linear  time,  or  at  least  in  time  0(n  \ogn).  In  attempting  to  enhance  the 
efficiency  of  our  technique  to  achieve  this  goal,  three  major  subproblems 
need  to  be  addressed. 

I.  Is  it  possible  to  calculate  the  relative  center  on  a  cut  in  0(n)  time?  Several 
new  ideas,  such  as  those  in  [Aea],  may  be  applicable  here. 

II.  The  second  subproblem  is  that  each  application  of  RELCEN  to  a  cut  of 
the  polygon  P  has  to  process  the  entire  collection  of  the  vertices  of  P,  even 
when  the  center  of  P  is  already  known  to  lie  in  a  smaller  subpolygon  of  P. 
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We  do  not  see  how  to  exploit  the  information  provided  by  previous 
applications  of  RELCEN  to  reduce  the  complexity  of  subsequent  applications 
of  this  procedure. 

in.  Finally,  one  needs  to  improve  the  complexity  of  the  final  stage  of  our 
algorithm,  involving  calculation  of  the  smallest  spanning  circle  of  n  given 
circles.  Can  this  be  done  in  linear  time? 
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